﻿using innovations.ml.core.solvers;
using innovations.util;
using innovations.util.exts.mathdotnet;
using innovations.util.exts.msft;
using MathNet.Numerics.LinearAlgebra.Double;
using MathNet.Numerics.LinearAlgebra.Generic;

namespace innovations.ml.core.models
{
    public class LogisticalRegressionModel : Model, IModel
    {
        public LogisticalRegressionModel(Matrix<double> x, Vector<double> y, bool includesOnesColumn) : base(x, y, includesOnesColumn) { }
        public virtual void ComputeCost(bool includeGradient = false)
        {
            int r = 1;
            double m = Y.Count; // number of training examples
            Vector<double> sigmoidValue = Sigmoid.Compute((X).Multiply(Theta));
            J = (1 / m) * (-Y * (MoreMath.Log(sigmoidValue)) - ((new DenseVector(Y.Count, 1) - Y) * (MoreMath.Log(r.PointwiseSubtraction(sigmoidValue)))));
            J = J + ((Lambda / (2 * m)) * (Theta.SubVector(1, Theta.Count - 1).PointwisePow(2)).Sum());
            if (includeGradient)
                ComputeGradient();
        }

        protected virtual void ComputeGradient()
        {
            double m = Y.Count;
            Vector<double> temp = new DenseVector(X.ColumnCount);
            Vector<double> sigmoidValue = Sigmoid.Compute((X).Multiply(Theta));
            temp = (1 / m) * X.Transpose().Multiply(sigmoidValue.Subtract(Y));
            if (Lambda != 0)
            {
                // We don't regularize the first column in X, 
                // so we can temporarily set the first theta element to zero.  
                // This elimnates the regularization term for first column.
                Vector<double> tempTheta = Theta;
                tempTheta[0] = 0;
                temp = temp + tempTheta.Multiply(Lambda / m);
            }

            // dm.Theta = temp doesn't work with L_BFGS solver in ALGLIB.
            // Therefore, we have to iterate through eack temp element 
            // and assign to dm.Theta.
            if (Solver.GetType() == typeof(L_BFGS) || Solver.GetType() == typeof(CG))
                for (int i = 0; i < temp.Count; i++)
                    Gradient[i] = temp[i];
            else
                for (int i = 0; i < temp.Count; i++)
                    Theta[i] = temp[i];
        }

        /// <summary>
        /// Used by ALGLIB.  All parameters are ignored.
        /// </summary>
        /// <param name="theta">ignored</param>
        /// <param name="j">ignored</param>
        /// <param name="gradient">ignored</param>
        /// <param name="obj">ignored</param>
        public virtual void ComputeCostAndGradient(double[] theta, ref double j, double[] gradient, object obj)
        {
            Theta = new DenseVector(theta);
            ComputeCost(true);
            j = J;
            // dm.Gradient = gradient doesn't work with L_BFGS solver in ALGLIB.
            // Therefore, we have to iterate through eack gradient element 
            // and assign to dm.Gradient.
            for (int i = 0; i < Gradient.Count; i++)
                gradient[i] = Gradient[i];
        }
    }
}